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Hydrodynamical stability of thin accretion discs: transient 
growth of global axisymmetric perturbations* 
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ABSTRACT 

The purpose of this paper is to explore how accretion discs manifest the phenomenon 
of transient growth on a global scale. We investigate analytically the time response of 
a thin accretion disc to particular axisymmetric perturbations. To facilitate an analyt- 
ical treatment we replace the energy equation with a general polytropic assumption. 
The asymptotic expansion of Kluzniak & Kita (2000), which extended the method of 
Regev (1983) to a full steady polytropic disc (with n = 3/2), is further developed and 
implemented for both the steady (for any polytropic index) and time-dependent prob- 
lems. The spatial form and temporal behaviour of selected dynamical disturbances are 
studied in detail. We identify the perturbation space which leads to transient growth 
and provide analytical solutions which manifest this expected transient growth be- 
haviour. Three terms (physical causes) responsible for the appearance of transient 
growth are identified. Two depend explicitly on the viscosity while the third one is 
relevant also for inviscid discs. The main conclusion we draw is that the phenomenon 
of transient growth exists in discs on a global scale. 

Key words: hydrodynamics - accretion discs 



1 INTRODUCTION 

It is widely accepted that accretion discs form whenever a sufficiently cool gas, endowed with a significant amount of angular 
momentum, is gravitationally attracted by a relatively compact object. This situation is quite common in astrophysics and, 
thus, the study, both observational and theoretical, of accretion discs has been quite intensive (see e.g. Pringle, 1981; Frank, 
King & Raine, 2002; for an overview). 

When accretion discs were theoretically proposed (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974) it had been 
recognized at the outset that the angular momentum transport necessary for accretion cannot take place unless some kind 
of enhanced effective transport process is invoked (typical values of the microscopic viscosity coefficient are by far too small 
to explain the observationally inferred accretion rates). As is well known, fluid turbulence induces anomalous transport and 
therefore angular momentum transport with the help of turbulent eddy viscosity has been proposed to operate in accretion 
discs. A detailed theoretical understanding of turbulence is still lacking and, therefore, the effective viscosity in discs has usually 
been approached in a phenomenological way: by parameterizing the effective viscous coefficient on the basis of dimensional 
arguments. This simple approach has been exceptionally fruitful, giving rise to successful interpretations of observational 
results (see Lin & Papaloizou, 1996; Frank, King & Raine, 2002; for reviews). 

The 'common knowledge' that turbulent flows arise by a transition from laminar flows, which become linearly unstable at 
a particular critical Reynolds number Re c , has largely been based on a number of well-studied experimental and theoretical 
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flows. However by the early 1990's the consensus within the astrophysical community was that this assertion appeared not 
to be applicable to thin Keplerian accretion discs (which are relevant to many astrophysical applications, in which the gas 
in the disc cools efficiently), because no definite linear hydrodynamical instabilities had been identified in such basic flows. 
Moreover, according to the famous Rayleigh circulation criterion, Keplerian shear flows are expected to be linearly stable. The 
magneto-rotational instablity, which has been shown by Balbus & Hawley (1991) to operate in such flows (when the fluid is 
electrically conducting and for not too large initial magnetic fields) , provides a linear instability mechanism for accretion discs 
and has since become accepted as the driver for magnetohydrodynamical turbulence and the resulting angular momentum 
transport. It has since become the operating paradigm in the astrophysical community that purely hydrodynamical turbulence 
in Keplerian discs is altogether ruled out - this conclusion is based largely on the paucity of hydrodynamic turbulent activity 
in numerical simulations (see Balbus, 2003; for a review and references). 

The question whether hydrodynamical turbulence can occur in accretion discs has however recently been re-opened in 
the astrophysical literature, close to a decade after Balbus , Hawley & Stone (1996) appeared to close this matter once and 
for all. The reason for this renewed interest in hydrodynamical stability of accretion discs appears to be three-fold. On the 
one hand, some instances of astrophysical systems have emerged in which the efficiency of the magnetorotational instability 
may be questioned (see e.g. Menou, 2000; Sano et al., 2000; Fromang, Terquem & Balbus, 2002). Secondly, a number of 
authors have applied to the accretion disc problem some recent ideas, familiar in the fluid-dynamical community, about the 
onset of turbulence and sustained dynamical activity in shear flows via the mechanism of transient growth (see below). Lastly, 
linear instabilities of shearing flows appear where gravity plays a catalyzing role. These stratorotational instabilities (Dubrulle, 
2004), which are distinct from instabilities giving rise to buoyant convection, have been identified from linear stability analysis 
(Molemaker et al. 2001, Yavneh et al. 2001, Dubrulle et al. 2004 Shalybkov & Rudiger, 2005) and have been proposed to 
operate in Keplerian discs. Yavneh et al (2001) demonstrate that this effect leads to significant non-linear activity in small-gap 
limit simulations of Taylor-Couette flows in a uniform gravitational field. It seems also possible that this instability may be 
responsible for the appearance of long-lived vortices in the localized hydrodynamical simulations of 3D Keplerian flows of 
Barranco & Marcus (2004). Though these are interesting, necessary and important avenues of exploration, these matters are 
not the subject of our attention in this work. 

The concept of transient growth (TG for short) of linearly stable modes and its possible role in giving rise to a bypass 
(of linear instability) transition to turbulence has been applied extensively to laboratory shear flows, which appear to become 
turbulent at Reynolds numbers that are significantly lower than Re c (given by the usual modal approach of linear stability 
analysis). In particular, for some archetypical flows (plane Couette, pipe Poiseulle), the theoretical critical Reynolds number 
is infinite, that is, these flows are found to be linearly stable. Yet such flow become clearly turbulent in the laboratory. The 
central idea is based on the fact the relevant operator, arising in linear stability analysis of shear flows, is generally non-normal 
(not commuting with its adjoint). Thus the usual procedure of finding the eigenvalues by solving the appropriate boundary 
value problem, may miss some important features of the linear problem. Even if all the eigenvalues imply stability (i.e. have 
negative or zero real parts), some initial perturbations may exhibit transient growth, because the eigenvectors of a non-normal 
operator may be non-orthogonal. The correct way to uncover such behaviour must then rely on solving the corresponding 
initial value problem. To save space, we shall refer the reader to a recent book (Schmid & Henningson, 2001) and a review 
paper (Grossmann, 2000), in which the idea that TG may trigger a transition to turbulence in shear flows is explained in 
detail, together with a comprehensive list of references. The astrophysical papers which will be mentioned shortly contain 
also, in their Introduction sections, short reviews on the history of this idea and its application to astrophysics. 

Among the papers in the astrophysical literature Ioannaou & Kakouris (2001) were the first to apply these concepts 
to the stability of accretion discs. They studied the global behaviour of incompressible two-dimensional (in the disc plane) 
perturbations, found copious TG for optimal perturbations and suggested that a turbulent state may be maintained by 
continuous random forcing. Most of the subsequent works utilized a local shearing box approximation and found linear TG in 
two dimensions (Chagelishvili et al., 2003) and in even in three dimensions (Tevzadze et al. 2003, Yecko, 2004). It was found 
that the third dimension may limit the very large amplification factors, but if the vertical scale is not too large the growth is 
still substantial. Yecko (2004) analyzed the linear three-dimensional problem and found that strong rotation essentially 'two- 
dimensionalizes' the most amplified transient disturbances whose maximal magnitude, for Keplerian shear, scale as Re 2//3 . In 
two very recent contributions Mukhopadhyayi, Ashfordi & Narayan (2005) and Ashfordi, Mukhopadhyayi & Narayan (2005) 
conducted detailed studies of the linear problem in various conditions and supplied estimates of the Reynolds number needed 
to sustain TG induced turbulence in accretion discs. They found that values of at least 10 5 — 10 6 are needed to overcome the 
stabilising influence of Keplerian shear (i.e. Rayleigh stable) and three-dimensionality. 

Recently, Umurhan & Shaviv (2005) employed a global asymptotic approach (utilizing the natural small parameter of the 
problem e = H/R, that is, the thickness of the disc relative to its typical radial scale) in the study of the dynamical evolution 
of global axisymmetric disturbances in an inviscid 3D disc. Such an asymptotic approach was introduced for the first time 

1 see Appendix A of Umurhan & Regev (2004) for a systematic discussion of this approximation 
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to the study of thin viscous accretion discs by Regev (1983) in the context of accretion disc boundary layers. It was later 
further developed and used in a remarkable analytical work by Kluzniak & Kita (2000) to solve for the steady structure of a 
polytropic viscous (using the conventional a prescription) axisymmetric 3D disc. This study revealed the presense of a steady 
meridional flow pattern with backflows for values of the a less than some critical value. This result and feature was confirmed 
by Regev & Gitelman (2002), who abandoned the polytropic assumption and included an energy equation (employing the 
diffusion approximation in the treatment of the vertical radiative transport). Umurhan & Shaviv (2005) allowed for time 
dependence and found algebraic temporal growth of global axisymmetric adiabatic disturbances, in a polytropic inviscid disc. 

The purpose of the present paper is to pursue this global asymptotic approach to pursue the time-dependent response of 
a general viscous thin accretion (a) disc and to identify the terms (mechanisms) responsible for the assortment of resulting 
behaviour, including any transient growth. The main departure that this work undertakes from previous investigations of TG 
in accretion discs (except that of Ionannaou & Kakouris, 2001 and Umurhan & Shaviv, 2005) is that we treat the dynamics of 
a sizeable part of the disc rather than a very small section centered about the disc's midplane (the shearing-box, e.g., Goldreich 
& Lynden-Bell, 1965; Balbus, Hawley & Stone, 1996), for which the boundary conditions are, strictly speaking, not known. In 
contrast to Ionannou & Kakouris (2001) we investigate here 3D (albeit axisymmetrical) disturbances and, as said before, this 
study generalizes the work of Umurhan & Shaviv (2005) to viscous flows and, as such, reveals the full range of possibilities. 
As we said, the primary tool utilized here in order to facilitate an analytical treatment is the asymptotic expansion, where the 
dependent variables and governing equations are expanded in powers of some small quantity (here the measure of the disc's 
'thinness', e). Exposing the resulting mathematical system to a set of perturbations, in which the extreme geometry of the 
disc structure is taken into account (i.e. its 'thinness'), we aim at obtaining an analytically treatable problem. To achieve this 
goal we also assume a polytropic relation between the pressure and density. Since this assumption was found to make little 
substantive difference to the steady meridional flow solution (Regev & Gitelman, 2002, and cf. Kluzniak & Kita, 2000), we 
see its use as justified here as well. The standard a prescription for the viscosity will be used here and thus we actually will 
be dealing with flows whose effective Reynolds number is Re ~ 1/a. 

The advantage of the asymptotic method and approach and the simplistic polytropic assumption is obvious - the treatment 
will be analytical and the responsible physical effects leading to any interesting dynamics may be transparently traced. It is 
clear, however, that the present analytic analysis should be ultimately complemented with a detailed and uncompromising 
numerical solution and a proper treatment of energy generation and transfer. 

As we shall see, the inviscid algebraic growth found in Umurhan & Shaviv (2005) will be replaced in this general viscous 
case by transient growth. The TG maximal magnitude will be seen to depend on a (i.e. the Reynolds number) and on the 
nature of the initial disturbance. It is still an open question as to 'if and 'how' non-linearities exploit the transiently growing 
disturbances to drive sustained turbulence. Given that this work is a linear study, this matter is not the subject of our 
discussions. We would like to remark, however, that repeated (every few hundred rotations or so) perturbative events can 
do the job of keeping dynamical activity alive (e.g., Ioannaou & Kakouris, 2002) and such events can not be ruled out in a 
real astrophysical accretion disc (e.g. due the presence of the second star in binary systems, variability of the circumstellar 
environment in young stars, various dynamical events in the vicinity of an AGN etc.). 

This paper is organized in the following way. In Section 2 we state all our assumptions, derive the basic non-dimensional 
equations and introduce the asymptotic expansions for the dependent variables that are used in this work. Although we are 
bound to repeat here some previously published work, we feel that explaining the global asymptotic approach to thin accretion 
discs (here including also time dependence) in a self-contained way may be useful for better understanding of our work here 
as well as perhaps for future investigations. In Section 3 the steady solutions of the equations are given and discussed. These 
solutions are then used to study the dynamical evolution by adding (in a suitable order) a time dependent disturbance to 
them. The dynamics of the disturbances is the subject of Section 4, where we present the results and discuss their nature. 
A global principle as well as detailed dynamical behaviour and some limits are discussed in an effort to uncover the origin 
of our primary result - transient growth. The case of a general polytropic index and some technical details of the analytical 
asymptotic procedure needed for obtaining the solutions are given in the Appendices. 



2 ASSUMPTIONS, NON DIMENSIONAL EQUATIONS AND ASYMPTOTIC EXPANSIONS 

Our starting point are the general Navier-Stokes equations 

/4^+p(V-V)V = -VP + ph + V-a, (1) 
g+V-( P V) = 0, (2) 

where V is the three dimensional velocity vector, p is the density, P is the pressure and b is the body force per unit mass. 
The Cartesian components of the viscous stress tensor, a are given by 
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dxk dx 



in which r\ is the dynamic viscosity coefficient (we neglect bulk viscosity because the processes discussed here are sufficiently 
slow to take place under thermodynamic equilibrium). 

In the context of accretion discs it is natural to present the equations in cylindrical polar coordinates with r, z and <j> 
being the radial, vertical and azimuthal coordinates respectively. Also, the disc matter's self gravity is neglected and therefore 
the body force derives from the gravitational potential of a central accreting object, whose mass is M, say. Thus we have 
b = — V<& with 

Regarding the viscosity of the accretion disc flow, we shall assume that the situation is such that the viscosity coefficient 
is greatly enhanced relatively to the microscopic one (see above in the Introduction). We will thus use one of the standard 
a-prescriptions (see below). If this effective viscosity enhancement results from ij actually being the eddy viscosity of an 
already turbulent flow, the variables of the fluid dynamical equations should be understood as mean quantities obtained by 
Reynolds averaging, a standard technique in treatments of turbulence (e.g. Monin & Yaglom, 1971). 

We proceed in a manner similar to that of Kluzniak & Kita (2000), hereafter KK, and Regev & Gitelman (2002), hereafter 
RG, by writing the equations in cylindrical coordinates (with the assumption of axisymmetry) in their non-dimensional form, 
which allows analytical asymptotic treatment. The difference is that here we shall not assume a steady flow: we retain the 
time derivative terms. This will enable us to analyze the dynamical evolution of deviations from the steady KK solution. 



2.1 The scalings and additional assumptions 

We consider here cold discs. Physically this means that the characteristic disc height (as measured from the midplane) is 
much smaller than its characteristic radius as measured from the accreting star's center. This feature is readily brought to 
the fore when we scale the dependent variables of the problem by their characteristic values (see Regev, 1983, hereafter R), 
which we shall denote by the 'tilde' sign. With R the characteristic radial scale of the disc (which is also a natural unit of the 
radial coordinate), the natural scaling for the angular velocity is the Keplerian value fl = ftk(R) = (GM/R 3 ) 1 ^ 2 . Also, the 
density will be scaled by a characteristic value (say that at the radius R and at the midplane of the disc): p. 

We assume throughout this work that the disc equation of state is polytropic - both in steady and dynamical states. This 
means that we can assume that the pressure and density always obey the relationship P = P(p) = Kp( 1+1 / n \ where n is 
the (constant) polytropic index and K is a constant. It follows that the typical scale of the pressure is P — P(p). Also, we 
choose the typical sound speed to be c s = \J P/p. The assumption that the disc is cold means that c s <C RO. which, in turn, 
allows for the definition of the vertical scale height of the disc, H = c s &- This is thus the natural 2 scale for the coordinate z. 
Formally speaking, we may define the small expansion parameter e to be, 

es ^ = ^« 1 . (5) 
RQ R 

This disparity of scales is exploited in what follows. The radial (V r ) and vertical (V z ) velocities are scaled by the typical sound 
speed c s and the angular velocity (fi = V^/r) by fi. Furthermore, we find that there is one natural choice for the temporal 
scaling, since a typical rotation time QT 1 and the vertical sound crossing time, H/c s , are defined to be identical. 

/.From here on all variables are assumed to be non-dimensionalized according to what has just been described. Thus the 
Keplerian angular velocity and the polytropic relations for the pressure and sound speeds are 

-=?'""■ <-^H i+ ;)H i+ ;)<-""■ (6 » 

To be consistent with previous works (e.g. KK) we make use of the non-dimensional sound speed c s as the dependent variable 
instead of the pressure P. As such it readily follows that 

ldP__ dcj 9c 2 s 
p dr dr ' p dz dz ' 

The standard a model of Shakura & Sunyaev (1973) is based on the assumption that the only non vanishing viscous 
stress component is a rt j, and it is proportional to the pressure. Following KK, we adopt this assumption and derive from it 
the form of viscosity coefficient, but include in the dynamical equations all the components of the stress tensor. In lowest 



2 Note that though a 'cold' disc implies H/R -C 1, a 'hot' disc, i.e. one where e ^> 1, does not imply that H/R 3> 1. Instead it just 
means that the structure starts looking more like a star with H/R ~ 1. 
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order in e the angular velocity of a disc is Keplerian and we get (with the dynamic viscosity coefficient scaled by c 3 H) the 
nondimensional relation (see also RG) 

2 oic 2 a 

where a is the viscosity parameter. 



2.2 Nondimensional equations and asymptotic expansions 



As stated in the Introduction, we seek axisymmetric solutions (in both the steady and dynamical cases) of the equations of 
motion. Denoting the non-dimensionalized radial and vertical velocities as u and v, respectively, the equations now appear as 
they do in KK and RG except that here we allow for time-dependence (the time unit is l/f2). Consequently time-derivatives 
are included and all the dependant variables are functions of r, z and t - 

3/2 



du 2 du du 2 

e— + e u— + ev- Or 

at or oz 



1 



e d / du\ 
p dz \ dz ) 

dv v 



dfl pu d(r 2 Q) dQ 

at r z or oz 
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dt dr dz 
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dz 



= or P? 



n dcj ld_ / dv\ _ ld_ / 2 ^ 

5r p dz \ dr ) p dr \ 3 dz 
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du 



dr ) 
del . 4 1 b 
P 

2 d ( r) d(ru 
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pr ar v ar / 
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dt r dr 



+ 



9(pv) 
dz 



= 0, 



(9) 
(10) 

(11) 
(12) 



where the gravitational potential of the central star has been expanded only up to second order in e. 

We seek asymptotic solutions to 19I12H by writing all the dependent variables in a power series of the small parameter e, 
i.e. we formally write 



U(r,z,t) = 



(13) 



where U = (u, v, Q, c^) T . We remind the reader that c s and p are not independent (on account of the assumption of polytropes) 
and are related through @. 

It has been shown before (see R, KK, RG) that in steady state the lowest order nonzero components of the meridional 
velocities are U\ and V2 and, in addition, the choice Q\ = pi = (and thus also c s i = 0) can be consistently made. Guided by 
these results (see also below) we retain in the expansions for Q, c? s , p and v only even powers of e (in the case of v starting 
from V2), while for the u expansions only odd powers are kept, starting with Ui, 



Q(r,z,t) = Q (r,z)+e 2 [tt 2 (r, z) + Q.' 2 (r, z,i)] + e 4 [Q 4 + Q' 4 (r, z,t)] +■ 

u(r,z,t) = e [ui(r, z) + u'i(r, z, i)] + e 3 [m(r, z) + u' 3 (r, z, t)] H 

v(r,z,t) = e 2 [v 2 {r,z) + v' 2 (r,z,t)] + e 4 [v4,(r,z) + v' 4 (r, z,t)] ^ 



c s (r,z,t) 



c s0 (r,z) + e 2 



c S 2(r,z) + c s2 (r,z,t) 



+ e 4 [cs 4 (r, z) + c' s4 (r, z, t)] + ■ 



(14) 
(15) 
(16) 

(17) 

(18) 



P(r,z,t) = p (r,z)+e 2 [p2(r,z)+p'2(r,z,t)]+£ 4 [p 4 (r,z)+p 4 {r,z,t)]-\ 

All O (1) quantities are assumed to be time-independent (in accord with the steady state results of RG and KK). As is also 
apparent, we have assumed no meridional flow at O (1) since, as said above, it has been shown before that ito = vo = vi = 0. 
Time dependence has been introduced into the expansions at O (e 2 J for the functions f2, p, c 2 and v while it appears at O (e) 
for u. At all orders in which time dependence is introduced we split the solutions up into a sum of a steady solution and a 
dynamical one, denoted by a prime, and this is true for all high orders as well. Thus, from now and on all the terms of the 
dependent variables expansions that are dependent on time are denoted by a prime, while those without prime are just space 
dependent, that is, steady. The implicit assumption in this kind of splitting is that the time-dependent part is a perturbation 
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(not necessarily infinitesimal!) on the steady state. In any case, the expansions Ansatz 1141 - 1181 1 can always be considered as 
a particular choice of a perturbation kind. 

The steady part of the solutions up to second order will be identical to the solutions obtained in KK and RG (for n = 3/2). 
In this work we shall limit the analysis of the evolution of the time-dependent disturbances to terms of order e 2 or lower. 



3 STEADY STATE 



3.1 0(1) 



All that is different between the results given this section and the results of KK is that all expressions are derived here in 
terms of arbitrary polytropic index n while in KK the specific value of n = 3/2 is assumed throughout their study. The lowest 
order equations are straightforwardly solved to yield solutions quite well-known in the literature (see, e.g., Hoshi, 1977). These 



c s o(r,z) 



h\r) 



2nr 3 



1 

r 3/2 ' 



Po(r,z) 



c 2 so(r,z) 



(19) 



The height of the disc h(r), i.e. where the quantities c 2 Q and po go to zero, is explicitly determined when solving the next 
order steady state equations. For values of \z\ > h, cj , po and all other associated quantities remain zero. Additionally the 
relationships between the sound speed, pressure and density, as well as their vertical gradients are given by, 



Po 
Po 



1 dP _ n dc 2 s0 
po dz dz 



z 



(20) 



It useful to note that we shall use in the equations to all orders only the zeroth order value of the viscosity coefficient, that 
is, will not actually consider it as an exapandable function, but rather as a prescribed function which is based on the zeroth 
order of the pressure. Thus 



r)o = ^aP r 



3/2 



(21) 



The steady state at this order of e has obviously a vanishing meridional flow (see KK) . 



3.2 0(e 



The steady-state equations at the next non-trivial order of e are, 
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Id, , d , 
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n dcj 2 | 4 1 
dz 3 po dz 



( dih ) 
po dz V dz J 

d ^ 



dz ) + po 



\d_ 

r dr 



( dui\ 



2d_ ( rj djrm) 
3 dz \ r dr 



(22) 
(23) 
(24) 
(25) 



Inspection of these equations show that they contain as solutions meridional flow, something which we will formalize below. 

KK found analytically the solutions at this order (i.e. equations 1221 - 1251 for the polytropic index n — 3/2. We discuss in 
considerable detail the solutions for an arbitrary polytropic index n in Appendix^ Here we would like only to stress some 
important points. 

First, as KK demonstrated, there is a value of a below which there is a certain amount of backflow in the disc, that is, a 
stagnation radius exists, beyond which the radial velocity, near the disc midplane, is directed outwards. 

Second, also as KK demonstrated, the disc height as a function of r, h(r), depends on the mass accretion rate, which in 
the appropriate units turns out to be of order e, M = eMi, and the polytropic index n. The disc's height h(r) can be derived 
in an asymptotically valid way only for radii significantly larger than r*, where r* is the zero torque position, i.e. the point 
at which the angular velocity (having the Keplerian value at lowest order) must pass a maximum on the way to matching 
it to some inner boundary value (the angular velocity of the surface of the accreting object, which is sub- Keplerian) . The 
point, if treated carelessly, will produce pathologies like zero disc thickness, a zero density and divergent velocities. To find 
the structure near r«, some boundary layer treatment, using singular perturbation methods like the ones employed by, for 
example, Bertout & Regev (1995) (see also references therein) or fully numerical approaches (see Popham & Narayan, 1995; 
and references therein), must be employed. Thus in this work, as well as in KK and RG, we are dealing with the disc only for 
radii sufficiently larger than r». 
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Third, these solutions do not allow for a detailed outer boundary condition (only an integral condition for the constancy 
of the total mass influx in the disc is employed). Thus, the KK solution and our generalization thereof (both steady and 
dynamical) are valid in regions of the disc that are not too close to its edges in the radial direction. 

Finally, it is an interesting result that in the expression for the disc height for arbitrary polytropic index (see IA5t the 
disc flares for n > 3/2, namely h(r) ~ r m with m > 1. 

The form of h(r) given in KK is recovered by substituting n — 3/2 into the general expressions (see lA5t given in Appendix 
A and reads as follows, 

The height of the disc at the fiducial radius r — 1 is denoted by hi and is 

fci = (2A)i(l-^,)*. (27) 



4 DYNAMICAL EVOLUTION 

The expansion procedure outlined in the previous section for the steady state solutions is extended here to dynamical problems. 
In other words, we dynamically evolve all equations at successive orders of e subject to the appropriate boundary conditions. 
Regarding the radial boundary conditions, we avoid the problems posed by them by considering only the portion of the disc 
with an inner radius ri n and an outer radius r max , as in the steady-state case (see KK). 

In the following, we discuss the vertical boundary conditions which we adopt. Let us say from the outset that the choices 
we make seem to us to be physically reasonable and, yet, allow for some analytical tractability. In steady-state, the disc lies 
between z = ±h(r), in which h(r) is as given in the previous sections. In this work we this consider the 'surface' of the disc 
to consist of the last fluid parcel which, in a steady disc, is at z — h(r). Technically speaking, the surface of the dynamically 
evolving disc should be given by h(r) + O ( e2 )> where the correction term accounts for the evolving surface. However, to the 
order up to which equations are sought in this work, it suffices to enforce the boundary conditions at z = h(r). 3 In addition, 
the equatorial dynamical behaviour allows two kinds of dynamical considerations according to whether the disturbances in 
the disc are sinuous or varicose. Varicose modes have even symmetry about the z = axis in all the quantities except the 
vertical velocity which has odd symmetry there. Sinuous modes have opposite symmetry of varicose modes. 

Firstly, thus, in this work we consider only varicose modes (the steady state meridional flow has this symmetry as well). 
However, it should be noted that when mode couplings begin to play a role in the evolution of the system (which occur at 
orders of e higher than that explored in this work), sinuous and varicose modes dynamically influence each other via the 
non-linear advection terms in the equation for the vertical velocity, that is, the vd z v term in 11 U provides this cross-mode 
interaction. 

Secondly, we allow the disc surface to flutter about freely in a Lagrangian sense. This, in turn, means that we require 
that the pressure on the moving surface vanish for all times. This, then, leads to the condition that the Lagrangian pressure 
fluctuation vanishes on this surface i.e., 

lk +eU lfr +V irz) P = Q > at * = ±h M- ( 28 ) 

Given the polytropic relation ©, the above boundary condition may be simplified with the aid of the mass continuity equation 
J3, to 



/ du dv 
\ dr dz 



at 2 = ±h(r). (29) 



The above is a restatement of the condition that there should be no thermodynamic work performed at the moving boundary. 

Thirdly, we posit that the viscous stresses on the surface of the disc are zero. This means that the product of the viscosity 
coefficient and the velocities (or the velocity gradients) tend to zero on a boundary. Since this is a free surface and the viscosity 
prescription given tends to zero at zero density, we will impose the more restrictive free-surface condition, 

lim r/rn ■ Vf2 -> and rjn ■ Vu -> 0, (30) 

z — >h 

where n is the normal to the surface S of the disc. In the event where the surface of the disc is approximately aligned with 
the normal to the vertical direction this condition becomes, 



3 In other words, the O (e 2 ) correction term to the disc surface does not alter the resulting disc dynamics at the lowest order for which 
they are determined. 
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lim rjr^^- — > and lim rj^- — > 0. (31) 

z->/t dz z->h dz 

In other words, we allow the gradient of the u and rQ velocities to diverge as one approaches the disc height h but the rate 
thereof should be slower than the rate at which r\ tends to zero. In the solutions that are developed for these modes (see 
below), we find that this condition is always satisfied and there are never diverging gradients of either u or VL. 

Lastly, we impose that the total vertically integrated mass flux through an annular section of the disc be fixed and steady. 
In turn this means that the vertically integrated mass flux through an annular disc section of dynamically varying quantities 
ought to be identically zero, i.e. 



27rr / pu'idz = 0, (32) 

J 'X 

in which i4(r, z,t) are the dynamically varying portions of the radial velocity at each order i of the e expansion. 



4.1 O (e 2 ) dynamical equations 

With the solution expansions implemented into the equation set (191 12t we find at the lowest nontrivial order describing the 
dynamical evolution of the disturbances (denoted by the primed quantities) the following, 



dtu[ 
d t v'i 



on rV 1 9 ( du '^ 
2ll Q, 2 r H — n— — 

po oz \ oz 

u'^ d(r 2 Q ) _j_d_ 
r 2 dr po dz 



dn 2 

'^7 



dc 2 s ' 2 4 19 
oz 3 po oz 



1 d 



dtp 2 = --^-(rpou'i) 



d 



dv' 2 
oz 

(pov' 2 ) 



r po dr 



djA_ 

' dz 



2_i_d_ 

3 po dz 



V - 



1 d(rv! x ) 
r dr 



(33) 
(34) 
(35) 
(36) 



r dr dz 

The polytropic equation state tells us that the relationship between the second order density and sound speed dynamical 
disturbances is 



a 1 P2 

CsO • 

n po 



(37) 



This algebraic relation is necessary to complete the linear PDE set 1331 - 1361 1. 

It is important to note that the above four equations support two sorts of modes, namely a pair comprising of purely 
vertical disturbances and a second pair containing disturbances in all the velocity components. The nature of these modes is 
detailed below. Such a classification is possible because equations I)33I34|I dynamically decouple from I|35I36|I . i.e. the quantities 
u[ and Q' 2 are not dynamically influenced by v' 2 and p' 2 , while the converse is not true. That is, disturbances in u[ and 
will drive a dynamical response in v' 2 and p 2 . This mathematical structure is reminscent of the behaviour of the steady-state 
equations (see section l3~2l and KK). 

To be more specific I|33I341 and 135I36|I may be combined to give two second order equations for u\ and r'>. 



= 0, 
Cv' 2 = [dtf+G]^, 

in which the operators are, 



(38) 
(39) 



V 

C 
T 

g 



„ i d d 

dt TT^TT 

po dz dz 



dl 



2 
c s0 



dz 2 



(n+1 



dc 2 s0 d „d 2 c 2 s0 



dz dz 



dz 2 



4 19 d n 
3 po dz dz 



2 19 9 19 9 

3 rpo dz dr rpo dr dz 
9 c 2 9 

dz rpo dr 

The zero Lagrangian pressure condition gives, at this order, 



du[ dv' 2 
dr dz 



0, at z = h(r), 



(40) 

(41) 
(42) 
(43) 

(44) 



which, we must reiterate, is at the unperturbed disc boundary. Also, the stress conditions are 
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r/r^ = and r,^- = 0, at z = htr). (45) 
az oz 

Finally, the integrated mass flux condition is equivalent to requiring 

27rr [ pou^dz = 0. (46) 

J-h 

We turn now to the presentation of analytic solutions of the system (1381391 . We shall first discuss, in section T4.2I the 
case u'i = (and thus also fl' 2 —0), i.e. solving only the homogeneous part of 1391 . The solutions of this equation constitute 
the above mentioned first pair of modes (purely vertical disturbances) and we shall refer to them as vertical acoustics. The 
solutions to the full set (1381391 in the general case will be discussed in section l4~3l This pair of modes will be referred to as 
driven general acoustics. The transient growth behaviour found in our solutions will be detailed in sect ion l4~4l A full exposition 
of the calculations leading to the solutions quoted and discussed in sections 14.214.41 is presented in in Appendices IB HB"2l In 
the next three subsections only the salient features will be given. 



4.2 Results: (i) Vertical acoustics 

The general solution to the inhomogeneous equation l|39|l can be written, in the form 

«2 = «h + v' v (47) 

where the index h stands for a solution of the homogenous equation and the index p stands for particular solution of the 
inhomogenous equation. As we have remarked earlier, the interesting property of the equations, which is a result of the disc 
being geometrically thin, is that the velocity and density fluctuations do not induce neither radial nor azimuthal motions 
at these lowest orders. This implies that homogeneous solutions of (1391 . (i.e. with the equation's RHS set to zero, because 
u'i = Q' 2 = say, see above) are perfectly acceptable. 
These solutions have the form (see Appendix IB11 

v' h = v' h (z,r)exp(vT) + c.c, with T=-^, (48) 

where the spatial eigenfunctions v' h are composed of the associated Legendre functions. This form of v' h already indicates that 
these solutions are inseparable in r and t. The temporal eigenvalues, which appear in complex conjugate pairs, are functions of 
a, the polytropic index n and the acoustic overtone (labelled by the integer index k, say). Thus the eigenvalues of the vertical 
acoustic modes should actually be written as a^{a, n). For instance, the eigenfrequency of the fundamental (i.e., k = 0) mode, 
is given by 



4 

a k= o(a, n) = of = —~ct±i 



16 2 2n + 1 



/16 



1/2 

, < a < 1 ; n > 1. (49) 



This shows that the fundamental mode comprises of decaying oscillations and all overtones show such a decay as well (see 
Appendix IB11 . In Figure we display the spatial structure of the eigenfunctions v' h and their vertical gradient dv' h /dz for the 
first three mode indices and for three polytropic indices. The structure of the eigenfunction is found not to be very sensitive 
to the polytropic index. 



4.3 Results: (ii) Driven general acoustics 

We turn now to the general solution to (1381391 by first focusing on 1381 . This equation admits an infinite set of eigenmode 
solutions in a way similar in character to the homogeneous solutions discussed in the last section. The general derivation and 
structure of these eigenmodes is detailed in Appendix IB2I For the sake of clarity, the discussion in the main text here shall 
focus mainly upon the dynamics associated with the fundamental mode of the operator V and we shall restrict our attention 
to the case n = 3/2. 

The solution to the dynamical radial velocity profile in an n = 3/2 polytrope is given by, 
ui = u(z, r)e CTT+¥; + c.c, with T = (50) 
in which the spatial eigenfunction w'(z,r) has the particularly simple structure 

«=Mr)^-l). (51) 

A(r) is an amplitude whose radial functional form is technically arbitrary and would be set by the initial disturbance condition. 
An arbitrary phase factor of tp is also introduced. The eigenvalues, denoted in this case by o, come in complex conjugate pairs 
given by 
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z/h z/h 

Figure 1. Structure functions for the homogeneous acoustic modes for a few values of the polytropic index and overtone number for 
values of a = 0.1. Plotted are the vertical mass flux M z = pov' h and Q z = Podv' h /dz. All Q z are scaled to their values at z = 0. 



a = ±i - fa. (52) 

Without loss of generality we need only consider the '+' solution since the phase factor tp takes care of the '— ' solution. The 
temporal behaviour of these modes is again one of decaying oscillations, but with frequencies given by the local rotation rate 
of the disc. Note that according to the prescription for h given in l|26[l . for values of r larger enough than the zero-torque 
radius r«, the functional dependence of the disc height with respect to r is well approximated by, 

h(r) ~ hir, ^~ ~ hi, 
dr 

where hi is determined by the mass flux rate and the value of a according to $27} or, for more general values of n, by l|A40 
and 1A6II . Since hi simply adds a multiplicative factor to all the dynamical quantities, it plays no role in the quality of the 
ensuing evolution and, as such, we set hi = 1 without any loss of general flavour. 

With this solution for u'i we may find a particular solution to the vertical velocity, v' p by solving 13911 directly. The details 
of this procedure are presented in Appendix [B2] but we highlight the major steps here. 

We posit the following form for v' p and p' p , 

v'p = v P (C,r)e aT + V p ((,r)(Te« T ) + c.c, (53) 
Pp = p P (Cr)e aT + R p ((,r)(Te° T ) + c.c, (54) 

where £ = z/h(r). Recall that we solve only for the fundamental mode of the system (again see Appendix IB2t and that 
there is obviously a dependence on a. It is important to notice explicit appearance of a multiplicative T term in the above 
expression. This temporal functional dependence is necessary in order to balance terms arising on the RHS of £23. In general, 
the form of these lowest order structure eigenfunctions is 

namely, it is polynomial in ( of the third degree with only odd powers in (. The four coefficients are functions of r and a, i.e. 
at — ai(r, a) and bi = bi(r, a) for i = 1,3. 

The solution to the density perturbation, p' p is straightforwardly obtained using l|51|l . 1530 along with po defined by l|19|l 
and then integrating 136H with respect to time. This gives 



Pp 



(i - O 



(56) 
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where, again, the coefficients are functions of r and a, that is, Cj = Ci(r, a) and di — dj(r, a) for i = 0, 2, 4. 

The polynomials appearing in the square brackets of the expressions for p p and Rp have only even powers of £. The 
detailed forms of the coefficients depend on the form of the generalized initial perturbation radial structure A(r) as well as 
on n. We avoid presenting them here because they are very long and cumbersome; and the details of the coefficients have 
no effect on the resulting TG since they describe only a particular type of disturbance. For the sake of simplicity all of the 
analytic results which will be presented below assume A = (that is, one particular form of the initial conditions of the 

radial disturbance) and n — 3/2. 



4.4 Results: (iii) Transient growth 

As is evident by inspection of 1531 . v' p and p' p exhibit TG. To be able to present it more clearly we consider the integral 
quantities 

£ a (r,t;a) = h(r) J (^p v' p 2 + l^^J d C E a {t; a) = £ a 27vr dr. (57) 

The dynamical quantity £ a (r,t;a) is to be interpreted as the acoustic energy (per unit area of the disc) consisting of the 
kinetic energy in the vertical velocity disturbances and the compression energy due to the density disturbances (see Section 
[SJ. The contribution of the (purely oscillatory) velocity components, resulting from the homogeneous part the disturbance 
equations, are left out when forming this energy integral. As £ a is a function of radial position r, one may form the quantity 
E a (t;a), i.e., the total disturbance acoustic energy (in this sense) of the disc, by integrating over the radial range in which 
these disturbances were assumed to exist, that is, between the inner and outer bounds r m i n (significantly above the inner 
zero-torque radius) and r max (significantly below the outer disc edge). 

In general, both £ a and E a are functions of the a parameter as well as the structure of the radial velocity perturbation 
A(r) (see equation I51H . As we have assumed (see Section [4.31 A(r) is taken here to be constant with respect to r, making 
the integrals defined in l|57|l analytically tractable. None the less the expressions, which have been verified with the aid of 
Mathematica 5.0, are very long and are not displayed here. Only their essential features, in particular their dependence upon 
a and r max , are addressed here. 

With A = e ,7r / 4 we find that £ a has the functional form 

-q2T 

£ a (r,t;a) = — g-^-T(T, cos 2T, sin 2T; a), (58) 

where T is a well-defined analytical albeit complicated function of its arguments. Thus, £ a depends on time only through the 
variable T = t/r 3 ^ 2 , which is the time measured in units of the local disc rotation period, at r, divided by 2-7T. Different values 
of r merely change the overall amplitude of the response as governed by the coefficient multiplying the function T. Otherwise, 
the evolution of £ a in the disc is self similar and is always decaying at long times T. Because of this we display below the 
evolution of £ a for different values of a at r = 1 only. 

Inspection of Figure [2] which plots against T the value of £ a (T) /£ a (0) for a range of a values, uncovers a very important 
property of TG. As a decreases the TG magnitude becomes more prominent (for example, it is up to ~ 3 orders of magnitude 
for a — 0.0025) and the maximum occurs at correspondingly later times. The time corresponding to the maximum amplitude 
Tmax is roughly ~ 5/ (8a) in this 4 case. The rise in energy is modulated by oscillations which arise from the fact that there 
are correlations between pairs of variables which contribute to the TG of £ a (T) (see Section 0. These variables oscillate in a 
frequency defined by the imaginary part of the eigenvalue (l~>2> and since products of pairs of variables terms depending on 
2T appear, so that the period is half the orbital period at the given radius. 

In Figure |3 we plot the behavior of E a for assorted values of a and r max . In all the cases displayed we have assumed 
that the inner radius, r m i n , is 0.05. For given r max , smaller values of a correspond to more dramatic instances of TG. The 
peak magnitudes are proportional to the inverse of a as is also the time when the peak occurs - a trend which is similar to 
that observed in the quantity £ a . It is very important to keep in mind that A(r) affects the spatial details but not the global 
time behaviour. Finally, we point out that the occurrence of the peak energy amplitude is delayed for larger values of r max , 
at given values of a, although the value of the energy at the peak time remains roughly the same. 

To get a sense for the expected qualitative spatio-temporal behavior of such transiently growing solutions we have 
computed the results for a case in which A(r) is not constant (and thus the possibility to evaluate the relevant integrals 
analytically can not be expected in general). Assume that the initial radial velocity amplitude has a Gaussian form centered 
around some fixed radius, 

(r-r ) 2 

A(r) = e m/4 e *f . (59) 
4 For the fc-th overtone, T max will be smaller by a factor of 1/fc 2 , see 1B12I 
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Comparative TG: Differing a at r = 1 



T = tr" 



Figure 2. The transient growth in the quantity £ a for the fundamental mode at r = 1, n = 3/2 and where A(r) = e 1 ""/ 4 , dA/dr = 0. The 
four curves correspond to a = 0.1, 0.01, 0.001, 0.0001 and are presented on a log-log plot. All curves are £ a scaled to its corresponding value 
at T = 0. The general rise time is proportional to the inverse of a as well as the maximum amplitude achieved. £ a is also characterized 
by oscillations with period 7r that sit atop the general transient trend. 

In our particular case we take the full width at half maximum given by A/ = 2.5 and the initial center of the disturbance 
ro = 3. Other parameters are n = 3/2 and a — 0.0025. As pointed out above, with this form for A we loose the ability to 
find an analytic solution and consequently resort to numerical evaluation. In Figure^Jwe show, in a contour plot in the disc 
meridional plane, the kinetic energy (per unit volume) contained in the vertical velocity disturbance v p , that is, poVp/2. The 
transient growth and decay of the disturbance is shown in the time sequence of figures displaying the spatial structure of the 
disturbance. We also depict how the disc surface moves in response to the imposed perturbation by solving the equation of 
motion for the boundary at this order. In the Figure time is quoted in terms of rotation times of the disc as measured at the 
radius r = 1. The two dramatic things to note is that, (a) there maximal rise in the disturbance energy is of almost three 
orders of magnitude and (b) the physical amplitude of fluctuations in the disc becomes quite large. Associated with the last 
effect is the emergence, as time moves on (by t ~ 100), of highly crenellated patterns on the disc surface. This crenellation 
is a result of the epicyclic frequencies causing the amplitude profile to wind-up on itself. Eventually, the winding is so strong 
that viscous effects take over and wipe the perturbations away by a few thousand rotation times. 

4.5 The ex — > limit: Algebraic growth 

The limit of zero viscosity, i.e. a — * 0, results in algebraic growth. We find that for n = 3/2, an arbitrary initial amplitude 
A(r) and the unapproximated disc height h(r), the vertical velocity fluctuation v' p simplifies to 



The density follows in similar suit. Notice that though the viscous effects vanish, the oscillation frequency associated with 
the local epicyclic vibrations remain. This solution appears as algebraically growing oscillations. This result is consistent with 
the more general (than that with the polytropic assumption) analysis made for discs with the fluctuations being adiabatic 
(Umurhan & Shaviv, 2005). 

It is also important to note that the limiting process of a — > is one that cannot be taken without caution and care. 
This is because this is a singular limit - substituting simply a — gives rise to a singularity in the expression for the steady 
state height h(r). Inspection of the solution to h(r) as appearing in IA5t shows that h diverges as a — » 0. However, one could 
formally take this limit holding the ratio M/a fixed (cf. IA5l and the definition of A in IA4t . namely 



lim v' = e 

o->0 





(60) 



< lim — < co, 

a— >o a 



(61) 
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Figure 3. Comparative TG behavior of the 'total disturbance energy' (E a ) and 'scaled disturbance energy' (Ea = E a (t)/E a (0)). In 
all plots n = 3/2, r m i n = 0.05 and A(r) = e 1 ""/ 4 Panel (a) depicts E a for different values of r max with a = 0.001. Larger values of r ma x 
simply correspond to longer lived phenomena before decay ultimately sets in at larger values of the disc radius meaning that the total 



energies are larger. In panel (b) Ea is shown 



E^ 



achieves the same overall maximum amplitude before decay. Similarly panels (c) 



and (d) depict E a and E a > for a = 0.1, 0.01, 0.001, 0.0001 at Tmax — 10. The trend for the E a to grow with decreasing a is evident, 
(d) we see that this growth, in terms of scaled energies, can be dramatic - by nearly three orders of magnitude. 



In 



Figure 4. Transient growth in the disturbance kinetic energy density in a disc where n = 3/2, a = 0.0025 and A(r) = e in / 4 e A / , 
with, in this example, ro = 3 and the full width-half-maximum Af = 2.5. Time is quoted in terms of disc rotation times at the radius 
r = 1 so that T = t. The magnitude of of the disturbance energy density is indicated by the various colours (whose meaning is indicated 
on the right of each panel) on a cut through the disc's meridional plane. A significant TG in the energy is observed (almost three orders 
of magnitude). Strong crenellation patterns appear in the fluctuating boundary of the disc before ultimate viscous decay takes over and 
wipes the pattern out. 



because the vanishing of the viscosity coefficient gives rise to a vanishing mass transfer rate. 



5 SUMMARY AND DISCUSSION 

This paper is a contribution to the ongoing investigation of the purely hydrodynamic response of accretion discs to small 
perturbations. Almost all of the recent studies, in which the phenomenon of transient growth has been found, were based on 
local (using the 'shearing box' approximation) analysis (see references cited in Section 0. While linear studies of this sort 
can be usually done analytically, pursuing the perturbation development into the non-linear regime calls for a fully numerical 
treatment. It appears that the existing numerical studies have not, as yet, achieved sufficient spatial resolution (or were 
not fully three-dimensional) to be able to definitely determine what is the ultimate dynamical fate of the transiently grown 
(sometimes by several orders of magnitude) disturbances (see, Yecko, 2004; Umurhan & Regev, 2004; Ashfordi et al., 2005; 
Mukhopadhyay et al., 2005). 

In an endeavor to investigate disturbances which are not limited to the shearing box (on which, typically, periodic 
boundary conditions are imposed), but still be able to treat the problem analytically or semi-analytically, we have resorted 
to asymptotic methods. Adopting such a method, previously applied (in R and KK) to steady thin (i.e. e< 1) axisymmetric 
viscous discs and extending it to include time-dependence, we are able to show that transiently growing disturbances naturally 
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emerge on the global scale as well. Moreover, we have found analytical expressions that help in understanding the underlying 
physics giving rise to the transient growth (see below). In achieving this goal a number of simplifying assumptions have been 
made. The unperturbed disc and the disturbances were assumed to be polytropic and axisymmetric and they were treated 
in a ring whose inner (outer) limit was significantly distant from the inner (outer) physical disc boundary. In addition, we 
chose to use the simplest initial conditions so as to excite the most fundamental modes, which we could readily analyze. As 
pointed out before these assumptions have largely been made to allow for analytical treatment in the framework of which 
we have derived a set of dynamic equations for successive orders of e. We stopped at the e 2 order, but the expansion can be 
carried out to higher order terms. The asymptotic expansion so defined brings out the salient features which are dictated by 
a geometrically thin cylindrical structure. Indeed, we were able to identify two kinds of modes and obtained the analytical 
structure of some simple initial disturbances. To O (e 2 ) we have determined an analytical solution to the dynamic problem 
and discovered that this set of equations manifests strong TG, irrespective of the polytropic index n. As we have noted, the 
severity of the TG is inversely proportional to a so that as the viscosity gets small, the TG becomes more dramatic. It seems 
then that TG is an intrinsic property of global axisymmetric disturbances of polytropic discs and, given the results of previous 
works done in the local 'shearing box' approximation as well as the two-dimensional (in the disc plane) global case (Ioannaou 
& Kakouris, 2002), we may reasonably expect it to be a property of all accretion discs. 

In spite of the simplifying assumptions we have implemented, and whose relaxation must be addressed in the future 
(see below), it is important to remark that the asymptotic expansion is not equivalent to the customary (in linear stability 
analysis) linearization procedure of infinitesimal perturbations imposed on a global steady state. Although, to the orders to 
which solutions have been determined the governing equations are linear, this procedure is intrinsically a finite-amplitude 
expansion (see Umurhan & Shaviv, 2005). Indeed, technically speaking, at each e" order, the equations are linear with 
inhomogeneous terms that result from the e n_1 order. Yet the aggregate solution, and the region of its validity, put it in the 
finite- amplitude realm. 

We turn now to the investigation of the physical sources for the TG that we have found in what we have called 'general 
driven acoustic modes'. To this end we consider the following energy density function 

2 

e a {r,z,t) = 1 p v 2 +2— P2- (62) 

The first term in this expression represents the kinetic energy contained in the time dependent part of vertical motion while 
the second term represents a sort of thermal energy in the form of an acoustic disturbance (Rayleigh, 1877). We have already 
considered integrals of a similar function, but only for the particular solution v' p (see Section^lJ. Here we shall be interested 
in the temporal behaviour of the volume integral of e a over the entire domain of interest r- m < r < r out ; \z\ < h. 

Using equations H35I36H to substitute for the relevant time-derivatives and converting the resulting integrals to surface 
integrals whenever possible in the usual way, we get 

I X = - X (p ^ o) dS+l L dS -"h &) 2 dS + Te ' (63) 

where Te represents volume integrals over the terms that are responsible for transient energy growth and its explicit form is 

The first term on the RHS of 1631 represents the mechanical work done on the disc by the bounding surface S. But, given 
that (a) the pressure of the moving boundary is zero and (b) the quantities in the integrand are odd functions of z, this 
surface integral vanishes. The second term represents the work performed on the surface by the viscous stresses associated 
with gradients of the vertical velocities. In a similar fashion the surface integral of this term also vanishes. Its general nature 
is not immediately obvious, but given the solutions we have obtained in this paper we always find that it vanishes, mainly 
because of the fact that r\ tends to zero (when approaching the disc surface) sufficiently fast. The third term is the volume 
integrated loss of energy due to viscous dissipation, i.e. the viscous losses. Since the integrand here is positive definite, this 
term is always a sink of energy and is the ultimate cause of decay in all the perturbations. Thus we concisely and transparently 
write the last relation as, 




e a dV — Te — [viscous losses]. 



Consider now the structure of the expression for Te in more detail. The first two terms in (1641 are explicitly viscosity 
dependent, while the third is not, although it may depend on a implicitly, via the a dependence of the disturbances. The 
first two terms consist of the cross-correlations between the vertical and horizontal viscous shear disturbances on one hand 
and the vertical velocity disturbance on the other hand. Regarding the third term, we notice that it follows from l|37fl that 
p'c 2 /po is proportional to the dynamical disturbance of the sound speed, i.e., c 2 2 ' (or, equivalently, of the temperature), so 
that the expression in the integrand may be rewritten as 



V d(ru' 2 ) 
r dr 



2 p djrpou'^ l dv 
po r dr J 



P2 c s 



(64) 



Hydrodynamical stability of thin accretion discs 15 



2 / 1 <9(rp wi) 



The volume integral of this expression thus depends on the correlation between the dynamical fluctuations in the radial mass 
flux and its radial derivative with the disturbances in the sound speed. Although we have chosen only disturbances which do 
not contribute at all to the total horizontal mass transfer rate (see the condition expressed in equation 1461 . the correlation 
integral consisting of the third term obviously need not vanish. In other words, the integral condition 1461 on the radial 
velocity dynamical disturbance in the radial velocity does not imply that the contribution to the work done on the fluid, as 
given in the third term, is zero. 

To assess the relative contribution of the first two viscous terms to above discussed third correlation integral, we calculated 
the relevant ratio and found that for the three values of a displayed in Figure [5] this ratio decreases monotonously with a. 
For example, for a — 0.001, it is less than 0.1% for all T > 0. Thus, for a — > the TG is totally dominated by the correlations 
embodied in the third term. This means that the algebraic growth found in the inviscid limit by Umurhan & Shaviv (2005) 
must be driven by this term. 

We conclude with some remarks on possible improvements to the asymptotic analysis of the sort done here and prospects 
for the future. Asymptotic expansions, when viable, are often very robust and provide a good approximation to the solution 
when truncation to only few leading terms is done. Obviously, when a term in the series becomes very large it may 'break its 
order', that is, become of the order of a previous term and as such make the expansion invalid in this region. In our expansions 
successive terms ratios are of O (e 2 ), and thus the procedure's validity should not severely be limited even up to a growth 
factor of ~ 100 or so (in the velocity or density perturbations). Still, the procedure, when carried to higher order, introduces 
corrections which are technically non-linear. Careful consideration must be undertaken in order to handle the response at 
these higher orders. This may entail treating the disturbance amplitudes for the lower order solutions (like A(r) in u[) as 
weakly non-linear governed by a second 'slow' time (e.g. the amplitude is instead written as A = A(r, r) where r = e 2 t) in a 
manner analogous to the treatment of non- linear thick polytropes (for example Balmforth & Spiegel, 1996). 

The approach used here may be generalized in a number of additional directions. Allowing for non-axisymmetric per- 
turbations, including the disc inner and outer boundary in some kind of boundary layer analysis and relaxing the polytropic 
assumption seem to be the most obvious generalizations. In such future studies the implications of the presence of backflows 
in the steady state should be examined as well. In any case, we have found the presence of TG in the simplest cases. It is 
difficult to imagine that it will be suppressed in the more general conditions although the effect of radiative energy losses on 
TG must be carefully examined. 

The question concerning the ultimate fate of the transiently grown perturbations and their ability to induce a sustained 
turbulent state in the disc remains open. In this context it is worthwhile to notice that since the TG decay times are of 
the order of hundreds of rotation periods, it is conceivable that accretion discs, which are usually not isolated systems, may 
experience recurrent external perturbations on such time scales and in this way the dynamical activity may be sustained. 

Extensive numerical calculations of as realistic as possible accretion discs are however needed to decide if TG per se 
may lead, through non-linear processes, to sustained turbulence. Such extremely high-resolution global three-dimensional 
calculations seem, however, to still be above the ability of the present computer power and it may be advantageous to also 
exploit sophisticated non-linear asymptotic methods to complement and guide them. 
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APPENDIX A: DISC FLOW SOLUTIONS FOR ARBITRARY POLYTROPIC INDEX 

We present here, for the sake of completeness, the steady-state solution, which can be found in KK only for a particular value 
of the polytropic index. 

Al Radial velocity solution and the disc height as function of r 

One begins by combining 1221231 into a single equation for Ui. After a little manipulation and upon using 121)1 we find: 



(h 2 - z 2 ) a 2 d_ 

2(n + 1) dz 2 Z dz 



. 1 I 2fi d ( 3 an 

+ 1 > Ui = 



r 2 po dr 



(»>£), (AD 



where h(r) is the (still unknown) disc height. The solution to 1A11 is 

, (~n- ah2 f (3n + l)[9(^ + l)(l-C 2 )+8(2n + 3)Q 2 ] 2r dh) 

{ ,U r 5 / 2 \ 9(n + l) 2 + 4(2n + 3) 2 q 2 h dr J 1 ' 

where £ = z/h. 

We are now in a position to determine the disc height as a function of the radius. First, we see that vertically integrating 
the continuity equation l|24|l . after multiplying it by 2nr, leads to a constraint whose physical meaning is that the mass flowing 
through the disc is independent on r. In particular, we are left with 

/ h ft 
dz— (rp oUl ) =0, (A3) 
■h 

because it is supposed a priori that the quantity po«2 indeed vanishes as z — > ±ft(r). The satisfaction of 1A3I is equivalent to 
saying that 

Mi = constant = — 2tt / dz(rpoUi), 
J -h 

where M\ is the mass accretion rate, which is one of the fundamental free parameters in this system 5 . 

The next step is the evaluation of the integral in the above equation using 1A2I . It leads to the following first order 
non-linear ODE for h(r) 

h 2n+3 , 2n + 3 r dh\ _ _ Mi T{n + 5/2) (2n + 2) r 



r 3 "+ 3 / 2 V 3n + l hdrj ^ a T(n + 1) 7r 3 / 2 (3n + l)' ^ 
Though the above is non-linear, its solution is remarkably easy to write down and reads 

h(r) = (2A)5^F3r^ ^/f^TK+s _ ( A5 ) 



5 Note that the minus sign in the above integral statement reflects the convention that a positive mass accretion rate Mi > implies 
physically that the integrated mass flux must be a negative in order for the material to flow inward. 
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This solution is shaped so that r, is the radius at which the disc experiences zero torque. At this radius our solution is singular 
(the thickness of the disc becomes zero and the radial flow infinite). For this reason we suppose that r» is always significantly 
smaller than the fiducial radius r = 1 and never consider these solutions near the zero torque point 6 Additionally, at the 
fiducial radius the height of the disc is given as 

hi = h(l) = (2A)5ST3 (i-^gsfcs . (A6) 

If hi is chosen to be some value, then there is a unique ratio of M /a which satisfies this constraint. 

Some things to note with regard to the disc thickness function: when n = 3/2 the height of the disc varies linearly with 
radius for r 2> r*; when n > 3/2 the structure of the disc is flared, meaning to say 4-— > as r gets large. The overall 
'thickness' of the thin disc, in which hi is its measure, is related to the accretion rate according to 

hi ~ Mf^ 3 . 

In the case where n = 3/2, the scale of the disc's thickness is proportional M^ 6 which is consistent with the results obtained 
by KK. 

In conclusion of this part of the Appendix we derive the general condition for the occurrence of backflow (see KK, RG). 
In other words, we find the highest value of a for which the radial velocity switches sign at some point in the midplane of the 
disc? This may be determined by setting to zero £ in IA2t and solving the result for a. After making use of JUJ we find that 
a and the position r at which the radial velocity in the midplane is zero are related through the equation 



There exists thus a critical value of a = that is, the one obtained for r — > oo in l|A7(l and it is given by 



3 \ln{n + 1) 



2 V (2n + 3) 2 ' 

Its significance is such that for values of a > there is no possibility of backflow in the disc. For values of < a < there 
are places in the disc, which includes the midplane, in which the flow is directed outwards. This flow transition, or stagnation 
point, bifurcates out of r = oo at the value of a = a x . For an n = 3/2 polytrope, ~ 0.685, which is the value predicted 
by KK. For values of a < the stagnation point, r stag occurs at 

r stag _81 (n + l) 2 (3n + l) 2 
r, 4 [9n(n + l)-2(2n + 3) 2 a 2 ] 2 ' V ' 

which, in other words, means to say that for a given value of a < tVoo, backflow occurs only for values of r > r s t ag - Finally, 
there is also a minimum radius r™ a g, below which there is never backflow. This occurs in the limit where a — > and is given 
by, 

'3n + 1 



(■£!)'. m 

In other words it says that r sta g appears in the a — > limit at the value r™ a °. In the case n = 3/2, this absolute minimum 
stagnation point appears at (121/36)r„ (see KK). 



A2 &,2,V2 and p-2 

The general solution to SI2 comes from inserting Ui into 123H and solving. We find that, 



a, = 



I, 2 



4 r 7/2 



(i-rt|§-.|i- 



4(3n + l)(l-(2n + 3)C )a 
9(n + l) 2 +4(2n + 3) 2 a 2 



(A10) 



Notice that the correction ^2 to the rotational speeds is non-zero in the limit of a — > 0. This merely represents the fact that 
at this order the pressure contributes something to the overall equilibrium of the system and this contribution makes the flow 
slightly SM&-Keplerian. Notice that unlike flo, the solution for £1? shows that the angular velocity is no longer constant on 
cylinders since 8^2/ d£ 7^ in general. 

Using 12m for po, together with 1A4L one may straightforwardly, though cumbersomely, integrate 1241 to yield a solution 
for the vertical velocity. The result is 



° None the less, it is clear that the infinities appearing near r* indicate that these solutions breakdown in that vicinity and that a 
boundary-layer type analysis will be required to cope with the failure. By restricting ourselves to be sufficiently far away from this 
pathological point we may proceed with the subsequent analysis. 
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(n + l)(2n + l)(3n+l)(l-C 2 ) 
(2n + 3)(9(n + l) 2 + 4(2n + 3) 2 a 2 



(All) 



The expression for the density is quite lengthly and we only present it for the case n = 3/2. The integration constants in 
deriving p2 were chosen so that c 2 p2/po — > as C, — > 0, 



P2 



(l-C 2 ) 3/2 h 5 
360^5r 13 / 2 



675(1 + C 2 ) + 16a 2 (460 + 528a 2 - 189C 2 
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r dh 
h dr 
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APPENDIX B: DETAILED CALCULATIONS FOR O (e 2 ) DYNAMICS 
Bl Vertical acoustics 



We begin this analysis with an exploration of the dynamics arising from the solution to the homogeneous part of equation 
0. (Bl) 



Recalling the arguments pointed out in the text, in this circumstance the vertical velocity and density fluctuations are 
completely decoupled from the horizontal (radial and azimuthal) motions. It is in this context that we consider the simplest 
vertical motions possible, Given that there are no horizontal velocities induced by these motions, the boundary condition I29H 
reduces to 

dv' h 



Po 



dz 







Assuming the solution Ansatz 14811 we find that (IBH reduces to 



1 / 8 ban 
^ n V + 9 n+ 1 



( 1 -C 2 )7^-2(n + l)C 7 --2n 



b' 2 + l 



•dC 



1 + 



9 n+1 



V h = 0, 



where for convenience we have introduced the height scaled variable £ = z/h. 

It is worth reviewing the behavior of possible solutions to the operator 1B3II near the singular points £ 
if we assume that the solution functionally behaves according to the form 



(B2) 



(B3) 



±1. In particular, 



v'h ~ (1 - 



C 2 )\ 



then an indicial equation analysis (Morse & Feshbach, 1953) shows that to leading order there are two values of A (corre- 
sponding to two linearly independent solutions) which are and — n. The former choice tells us that there is, in fact, a solution 
described by a regular series expansion near ±1. As for the latter solution, since it predicts that v' h diverges like (~ n the 
product Podv' h /dz none the less remains constant. However this solution must still be rejected because it would cause the 
zero Lagrangian pressure condition on the surface to be violated, i.e. (14411 would not be satisfied. 

As it turns out, the solution to 1B3II is composed of the associated Legendre polynomials P and Q (Abramowitz & Stegun, 
1972) as follows, 



with 
A(a,n,a) 



(i-C) 



i/2 



Mr) [P2(OQl(o)-Ql(OP2(o)] - 

2 8n(n + l)(l + a 2 ) 



-1 + 



1+ (2n + lY 



n + 1 + |nacr) 



1/2 



(B4) 



and this solution satisfies the condition 144H . The z independent constant A(r) is entirely arbitrary and is determined by the 
initial conditions. In order to satisfy the boundary condition at z = h a somewhat complicated expression emerges for the 
eigenvalue condition. However, it turns out that for special cases of the polytropic index the quantization condition for these 
modes may be easily written down, namely, for n a positive integer or half odd integer it is required that, 



A = 2k + n + 1, 



(B5) 



where k is any integer k > 0. In this sense k labels the acoustic overtone of the mode with k — being the fundamental. Given 
the definition of A above, for every k there are two values of a which are always decaying. In considering all fundamental 
modes, which we label with <r F , we find that, 
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for a < a Q 



(B6) 



It is interesting to note that for a < aw , given by, 



Q^osc — 



9 / 2n + 1 



the modes are complex conjugate pairs of decaying oscillations. For a > a oac the modes are strictly exponentially decaying. 
For example, for n — 3/2, a osc — (9/4) ■ -y/8/3 ~ 4. Furthermore, in the zero viscosity limit, the oscillation frequency is given 
by 



a F (a — *• 0) = ±i 



2n + 1 



B2 General Driven Acoustics 

We now turn to solving the coupled equations I38I39H . Given the structure of these equations we must begin with developing 
a solution to 1381 and then using this in We suppose, as we did in the previous section, that the solution to 1381 . which 
is a homogeneous equation, has the following form, 

u'i = u(z, r)e aT + c.c. (B7) 

We write the eigenvalue as a in order to distinguish it from the 'vertical acoustics' eigenvalue, a, discussed in the previous 
section. After a little manipulation we find that 13811 turns into the far simpler ODE 



+ 1 }u = 0, (B8) 



2(n+l) v W 9C 2 d( 

Before presenting the Ansatz for the solution it is worthwhile to consider the array of solution behavior of l|B8ft near the 
singular points f = ±1. Assuming that near f = 1 the behavior is u ~ (1 — C, 2 )^ , an indicial analysis (Morse & Feshbach, 
1953) reveals that there are four choices for j3, namely, 0, 1, —n, 1 — n. The solutions with /3 = —n, 1 — n are rejected on grounds 
that they violate the boundary condition at £ = ±1 (they essentially diverge). In other words, since r\ is proportional to Po, 
then as one approaches the boundary, rj(du' /dz), rjr(dfl' /dr) 7^ there for these two values of j3. It would then mean that 
the stress conditions I45H would be violated there and so these solutions are ignored. The remaining two possibilities suggest 
that a regular truncated series expansion is a satisfactory means of expressing the solution to u' . 

Thus, for varicose modes we assume the Ansatz 

«' = fi W = E^MC 2(fc+1) - 2m , (B9) 

where k labels the eigenmode and is an integer greater than or equal to zero. Correspondingly, the eigenvalue, a, would be 
labelled by k too. Aside from a generalized characterization of the eigenvalues for arbitrary k (see below), from here on we 
restrict ourselves to considering only the fundamental mode, i.e. k = and thus, 

6' = fi {0 ) = A(rK 2 + B(r), (BIO) 

where the constants have been redefined for simplicity. Generally speaking, insertion of the solution form 1B91 into IB8t 
requires one to recursively solve for the values of the coefficients Ak in the usual procedure called for when developing series 
solutions to ODE's. The first of these equations is always a relationship for the eigenvalue o. For example inserting IBlOt into 
1B8I I gives 



A 



( 22n + 3\ 2 , a2 ( n „2 a \ .2 a ( 2 2n + 3\ n „ 

(o + a +1 C + o\oB-A — A 1(7 + -a + B = 

V 3n+l/ s V 3n+17 3n+lV 3 n+l J 



and by setting to zero the coefficients of each power of £ one gets 

2ti + 3 B(r) 1 . . 

a — a± — ±i — a , , , ' = . (Bill 

n + l' A{r) 2n + 3 K ' 

In other words, there are the two possible temporal responses (denoted by <r±), to this type of disturbance, however, the 
coefficients are the same for both. We note that the eigenvalue form indicates that these modes exhibit decaying oscillations. 
For general values of k, the frequency response obtained by inserting 1B9I into 1B8II . is given by 

(fc + l)(2n + 3 + 2fc)2 ,„ , 

a k =a(k) = ±i- ^ j—TTi " o Q ' B12 

(n + l) 3 
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showing that these modes become increasingly damped in proportion to k 2 in the large k limit. The general solution for an 
arbitrary disturbance can thus be written as 

oo 

u'i = y^a (fc) (C,r-)e' TfcT +c.c. 

fc=0 

We note also that the integrated radial flux of mass for the perturbation vanishes, namely that, 

poru'dC = 0. (Bf3) 



-i 

and this is generally true for any integer mode k. The consequence being, as will be seen below, that the radial structure 
of the perturbation amplitude, namely the quantity A(r), remains a freely chosen function, usually guided by some initial 
condition. 

Now we turn to expanding the particular solution to 139^ armed with the solutions ijTTFt . llBTnjl and ijBTTjl . Inspection of 
the operator structures of C, T and Q suggests that the solution to the particular vertical velocity, v' p , has the form 

oo 

^ = I](«(fc)(C,Oe CTT + V (fe) (C,r)r e CTT ) +c.c, (B14) 
in which 

k + l fe + 1 

E„ / \ A2(fe + l)-2m + l t> V""* r\ I \/-2(k+l)-2m+l / D ,r\ 

£ifc,2m(r)C k ' , V (k) = Sfc, 2m (r)C k ' ■ (B15) 

m=0 m — 

We note that these solutions have terms which manifest TG accompanied by oscillations due to the form te~ at e 1 ^ 1 type of 
behavior. Recalling that we are here only considering the fundamental mode, the above solution form is written as 

v {0) = HooMC 3 + H 02 (r)C, V (0) = 6oo(r)C 3 + 6 02 (r)C. (B16) 

Using l|B16|l together with the solutions l|B7|l . ljB10fl in l|39|l and after some lengthly algebra (which was verified using Mathe- 
matica 5.0), the coefficients straightforwardly follow. Unfortunately, their expressions are very long and we have omitted their 
presentation here. This solution for v' p along with the solution for u' ± , when used in 1361 immediately yields the solution to 
p'i- Its expansion is similar to that of of v' p , namely, we say that p' 2 ~ p' p + p' h , representing its particular and homogeneous 
contributions. Then it follows that 

p; = p'(C,r)e CTT + J R'(C,r)Te CTT +c.c, (B17) 
in which 

p' = P{k) = g H^ m M< 2(fc+1) - 2 "\ R' = R (k) = E e^ m (0C 2(fc+1) " 2m . (Bis) 

m— m— 

The superscript (p) appearing in the above expressions is meant as a label to distinguish these functions from the corresponding 
ones for the vertical velocity mode, v' p , cf. 1B15L The lowest density fluctuation eigenmode is 

p (0) = H^(r)C 4 + H^(r)C 2 + H^(r), R (0) = 6« (r)( 4 + 9^(r)C 2 + Q${r). (B19) 

The solutions quoted in the text have this form. We note also that the boundary conditions on the vertical velocity 
components, namely, 



2\ "+ 1 



c 2 ) 



1 dv' p 1 dru^ 



h(r) d( r dr 
as £ — ■> 1, is satisfied by these solutions 



0, (B20) 
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